#include "defs.h"

int main(int argc, char** argv)
{
	int n = 6;
	
	vdouble res (n-2);
	vdouble x (n-2);
	vdouble diag (n-2);
	vdouble diag_otra (n-1);
	
	forn(i, n-2)
	{
		diag[i] = 1;
		diag_otra[i] = 1;
		res[i] = i+1;
	}
	
	res[0] = 2;
	diag_otra[0] = 0;
	diag_otra[1] = 2;
	diag_otra[n-2] = 1;
	
	cout.precision(12);
	cout << "Matríz tridiagonal " << endl;
	
	forn(i, n-2)
	{
		forn(j, n-2)
		{
			if (i == j) 
				cout << diag[i];
			else
				if (j == i + 1 )
					cout << diag_otra[j];
				else
					if (j == i - 1)
						cout << diag_otra[i];
					else
						cout << "0";
						
			cout << "\t";
		}
		cout << endl;
	}
	
	cout << endl << "Res" << endl;
	
	forn(i, n-2)
		cout << res[i] << " ";
	cout << endl << endl;
	
	vdouble coef_c (n-2);
	vdouble coef_d (n-2);
	
	coef_c[0] = diag_otra[1]/diag[0];
	coef_d[0] = res[0]/diag[0];
	
	cout << "c'[0] = " << coef_c[0] << endl;
	cout << "d'[0] = " << coef_d[0] << endl;
	
	fornx(i, 1, n-2)
	{
		coef_c[i] = (double) diag_otra[i+1]/ (double) (diag[i] - coef_c[i-1] * diag_otra[i]);		
		coef_d[i] = (double) (res[i] - coef_d[i-1] * diag_otra[i])/ (double) (diag[i] - coef_c[i-1]*diag_otra[i]);
		
		cout << "c'[" << i << "] = " << coef_c[i] << endl;
		cout << "d'[" << i << "] = " << coef_d[i] << endl;
	}
	
	x[n-3] = coef_d[n-3];	
	fornd(i, n-3)
	{
		x[i] = coef_d[i] - coef_c[i] * x[i+1];		
	}	
	
	cout << endl;
	
	forn(i, n-2)
	{
		cout << "x[" << i << "] = " << x[i] << endl;
	}
	
	cout << endl;
	
	vdouble c (n);
	
	c[0] = 0;
	c[n-1] = 0;
	
	forn(i, n-2)
	{
		c[i+1] = x[i];
	}
	
	vdouble h (n-1);
	vdouble a (n-1);
	
	/* Llenar los vectores. */
	forn(i, n-1)
	{
		h[i] = 0.2;
		a[i] = 0.1;				
		cout << "h[" << i << "] = " << h[i] << endl;		
		cout << "a[" << i << "] = " << a[i] << endl;		
	}
	
	cout << endl;
	
	/* Obtener el resto de los coeficientes en función de los c_i. */	
	vdouble b (n-1);
	vdouble d (n-1);
	
	forn(i, n-1)
	{
		d[i] = (c[i+1] - c[i]) / (3.0 * h[i]);
		b[i] = a[i]/h[i] - h[i]/3.0 * (2*c[i] + c[i+1]);
	}
	
	forn(i, n-1)
	{
		cout << "c[" << i << "] = " << c[i] << endl;
		cout << "d[" << i << "] = " << d[i] << endl;
		cout << "b[" << i << "] = " << b[i] << endl;
	}
	
	return 0;
}
